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Dark-dark solitons and modulational instability in miscible, two-component 

Bose-Einstein condensates 
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We investigate the dynamics of two miscible superfluids experiencing fast counterflow in a narrow 
channel. The superfluids are formed by two distinguishable components of a trapped dilute-gas Bose- 
Einstein condensate (BEC). The onset of counterflow- induced modulational instability throughout 
the cloud is observed and shown to lead to the proliferation of dark-dark vector solitons. These 
solitons, which we observe for the first time in a BEC, do not exist in single-component systems, 
exhibit intriguing beating dynamics and can experience a transverse instability leading to vortex 
line structures. Experimental results and multi-dimensional numerical simulations are presented. 

PACS numbers: 03.75.Kk, 67.85.De, 47.40.x, 03.75.Lm, 05.45.Yv 
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Superfluids are a robust model system for the inves- 
tigation of nonlinear fluid flow. Governed by an un- 
derlying macroscopic wavefunction, superfluids can dis- 
play a large variety of nonlinear wave phenomena in the 
context of matterwaves. In Bose-Einstein condensates 
(BECs), nonlinear structures including solitons, vortices 
and vortex rings have been the focus of intense research 
efforts [l|, [2[. In this work, we investigate the regime 
of fast counterflow between two distinguishable superflu- 
ids in a narrow channel and observe dynamics leading 
to novel structures. Modulational instability (MI), in 
which small perturbations to a carrier wave, reinforced 
by nonlinearity, experience rapid growth [3[, plays a key 
role in the dynamics. In many nonlinear systems, MI 
leads to the breakup of periodic wavetrains, as in suffi- 
ciently deep water [4| , as well as the formation of localized 



structures in optics [5j and BECs [6J. In our case, MI- 
induced regular density modulations, formed throughout 
the BEC, lead to the emergence of a large number of 
beating dark- dark solitons. These solitons-which exhibit 
periodic energy exchange between the two condensate 
components [J] -are a generalization of static dark-dark 
solitons [3|- They are distinctly different from all pre- 
viously observed solitons in BECs, including dark-bright 
solitons which were generated in a two-component mix- 
ture by marginally critical counterflow-induced MI near 
a density edge [9|. We perform three-dimensional (3D) 
numerical simulations to corroborate this interpretation 
and furthermore identify a subsequent transverse insta- 
bility resulting in multi-dimensional structures such as 
vortex lines (see [10] for the scalar counterpart). 

We study superfluid counterflow with an experimen- 
tal system consisting of BECs with typically 8 x 10 5 
atoms of 87 Rb. The BECs are confined in a cigar shaped, 
far-detuned optical dipole trap with measured trap fre- 
quencies of 27rx{1.5, 140, 178} Hz with a horizontal 
weakly confined axis. By starting with all atoms in the 
\F,rriF) = |1, — 1) hyperfine state and transferring about 
50% of the atoms to the |2, — 2) state via a 1 ms long 



microwave sweep, a perfectly overlapped two-component 
mixture is created. The predicted scattering lengths for 
these states 11( imply that this mixture is miscible [12|, 
which is also supported by our experimental observation 
of no phase separation for an unperturbed mixture of 
these states. To induce relative motion between the com- 
ponents, an external magnetic gradient is applied along 
the elongated (axial) direction. The gradient pulls atoms 
in the |2, — 2) state to the left and those in the |1,— 1) 
state to the right. The atoms are imaged using a free 
expansion imaging procedure. Each experimental image 
shows an upper cloud consisting of the |2, — 2) atoms af- 
ter 7 ms of free expansion and a lower cloud consisting 
of the |1, — 1) atoms after 8 ms of free expansion. Both 
clouds are imaged during the same experimental run. 

Experimental data showcasing the formation of a very 
dense counterflow-induced MI pattern are presented in 
Fig. [TJ In the presence of a 10.4 mG/cm axial gradient, 
gradual pattern formation starts after 70 ms of smooth 
evolution (Fig. QJa, b)). We first observe pattern for- 
mation in non-central regions where the two condensates 
have differing densities (Fig.QJb)). This is due to the de- 
pendence of the critical velocities for counterflow-induced 
MI on the two component density ratio, being largest 
when the densities are equal [5|, |9[. After about 25 ms, a 
very dense and regular MI pattern fully develops, filling 
the entire BEC (Fig. 0Jc)). The modulations in the two 
components are offset in the axial direction in a staggered 
way such that one component fills the depressions in the 
other. Under the continued influence of the axial gra- 
dient, the regular pattern of Fig. []Jc) quickly becomes 
uneven and irregular. Alternatively, if the gradient is 
switched off after the MI pattern has fully developed, 
we frequently observe the formation of black dots such 
as those marked by the arrows in Fig. QJd) which might 
indicate the generation of vorticity. We note recent theo- 
retical work suggesting that counterflow-induced MI may 
be used to generate quantum turbulence |14j . 

Imparting slow counterflow conditions, implying slow 
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FIG. 1: Counterflow induced MI in the presence of a strong 
magnetic gradient of 10.4 mG/cm. Evolution times (a) 10 ms, 
(b) 70 ms, (c) 95 ms. (d) After MI onset, the axial gradient 
is turned off, followed by a trapped evolution time of 20 ms. 



MI onset in the quasi- uniform background, we previously 
generated a dark-bright soliton train emanating locally 
from a density edge [9j]. In contrast, the fast counter- 
flow considered here leads to rapid MI onset and pattern 
formation throughout the condensates. 

MI theory agrees quantitatively with the experimen- 
tally observed patterns as we now explain [Fig. [2]. For 
a uniform counterflow, the onset of MI corresponds to 
a complex sound speed (see [15]) and exhibits a pre- 
ferred wave number, /c max , corresponding to the maxi- 
mum growth rate ^ max , both depending on the counter- 
flow speed. Unfortunately, our imaging procedure does 
not allow us to determine the counterflow speeds exper- 
imentally. However, we can take two independent the- 
oretical approaches, described below, to extract the on- 
set velocities from our experimental data. The fact that 
these two independent approaches lead to consistent re- 
sults gives quantitative credence to the theory. First, we 
use the analytical theory in [5|, [9| to calculate the coun- 
terflow speed Vfit whose corresponding /c max equals the 
experimentally observed pattern periodicity at the trap 
center where the densities are assumed equal (solid, black 
curve in Fig. [2j). In a second, independent approach, we 
assume spatially uniform counterflow whereby the ap- 
plied gradient leads to unimpeded acceleration of each 
component (calculated from the atomic magnetic mo- 
ment and the magnitude of the applied gradient). Us- 
ing this simple model, experimentally determined onset 
times for MI are converted to relative speeds at the on- 
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FIG. 2: (Color online) Consistent predictions of counterflow 
speed based on wavelength and onset time measurements of 
MI as a function of applied gradient. For details see text. 



set of the MI pattern (dashed blue curve in Fig. [2]). The 
dash-dotted red curve in Fig. [2] uses the same, uniform 
counterflow model but shifts the measured MI onset time 
by — 1/^max- Subtracting this time accounts for the de- 
velopment of the instability and leads to a better ap- 
proximation of the true relative speed that sets the pat- 
tern periodicity. The resulting curve interpolates the two 
models. The lowest, dotted curve is the predicted critical 
speed in the condensate center (v cr = 0.16 mm/s) demon- 
strating fast counterflow. Despite the approximations 
made, the curves exhibit agreement for small to moder- 
ate gradients, suggesting that the observed dynamics are 
theoretically described by counterflow-induced ML Dis- 
crepancies at large gradients are likely due to the large 
accelerations involved and spatial nonuniformity. 

We now investigate the dynamics of the MI onset by 
using a smaller gradient of 1.4 mG/cm so that /c max is 
reduced relative to Fig. [TJ enabling better experimental 
observation of individual features (Fig. [3]). After smooth 
counterflow, MI sets in across the BEC leading to a reg- 
ular array of dark-dark solitons (Fig. [3fa, b, d, e)). In 
accordance with theory and our numerics (see below), 
the dark-dark solitons exhibit a dynamic beating as seen 
by comparing the integrated cross sections of Fig.[3jd,e), 
noting the order of the notch and bump feature in each 
component. While our destructive imaging technique 
does not allow us to determine the exact beat frequency, 
our 3D numerics indicate a timescale of fifteen millisec- 
onds per period [15] . The dark-dark solitons we observe 
here are new and distinct from the dark-bright solitons 
that have been observed previously in BECs [9|, LL6|, LL7j , 
being distinguished by their far field conditions and dy- 
namics. To facilitate a comparison, an example dark- 
bright soliton train, seeded at the condensate interfaces 
and generated by slow, marginally critical counterflow 
[9|, is shown in Fig. [3^c, f). A dark-bright soliton con- 
sists of a dark notch in one component, filled by a local- 
ized density bump of the second component. In contrast, 
the beating dark-dark soliton asymptotes to nonzero den- 
sities in both components and dynamically changes its 
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FIG. 3: (Color online) (a, b) Dark-dark solitons as a result 
of MI after applying a gradient of 1.4 mG/cm for 350 ms and 
380 ms, respectively, (c) Formation of dark-bright solitons 
when a small magnetic gradient of 0.2 mG/cm is applied. The 
gradient is left on for 1000 ms before the start of the expansion 
sequence, (d-f): Zoomed-in view of boxed region in (a-c), 
respectively. Red solid lines are integrated cross sections of 
|2, —2) state, black dashed lines of |1, —1) state. 



shape, with each component possessing a density bump 
adjacent to a notch which alternate their relative posi- 
tions in time [see also Fig. [6] below]. 

The dynamics are well reproduced by 3D numerical 
simulations [15j of the vector, mean-field Gross-Pitaevskii 
equation with initial conditions and parameters corre- 
sponding to the experiments in Figs. QJ a-c) and[3fa,b). 
As with experiment, a smooth, accelerating counter- 
flow develops due to the axial field gradient. Dark- 
bright solitons form at the edges of the condensates until 
the rapid growth of large scale modulations is observed 
(Fig. Sfa,b)). For moderate gradients in Figs. 3Ja,c,e), 
these modulations rapidly develop into a number of lo- 
calized, essentially one-dimensional (ID) beating dark- 
dark solitons with initial approximate spacing 27r/& max . 
Continued evolution results in interactions and eventual 
solitary wave transverse breakup at about t = 600 ms. 

For the strong gradient case, our numerics show the 
development of axial modulations by about t = 125 ms 
with an initial ID structure. In contrast to the moderate 
gradient regime, these structures rapidly undergo decay 
due to transverse modulations which leads to the forma- 
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FIG. 4: (Color online) Integrated densities from 3D numerical 
simulations. (a,c,e) correspond to Fig. [3ja,b) at t — 421 ms 
with zoom in (c) and line plot (e) of dark-dark solitons. (b,d,f) 
correspond to Fig. [TJa-c) at t = 133 ms with zoom in (d) 
and a phase plot along the vertical z = plane (f) showing 
two vortex lines with oppositely oriented 2ty phase winding. 
The vertical axes of (a-d) span 16.7 microns incorporating a 
vertical offset of 8 microns between the clouds. 



tion of columnar 2D vortex lines, Fig.|4jb,d), exhibiting a 
2tt phase winding around their core, Fig.H^f), and a uni- 
form structure along the direction of view. The numerics 
also show vortex lines oriented along the orthogonal, hor- 
izontal radial axis. In analogy to the scalar case [lOj, we 
interpret this behavior as a transverse instability that de- 
pends on the relative speeds of the two components, their 
densities, and the transverse confinement strength. 

Dark-dark solitons can also be observed in other set- 
tings, e.g. during the mixing of two initially phase sepa- 
rated components. An experimental result is showcased 
in Fig. [5] We start from the phase separated situation 
in Fig. Efa,c) which forms after initially overlapped com- 
ponents experience an axial gradient for 10 sec. When 
the axial gradient is suddenly switched off, the two com- 
ponents interpenetrate, first forming a smooth and ex- 
tended overlapped region. After some evolution time, in- 
dividual dark-dark solitons appear (Fig.[5fb, d)) exhibit- 
ing approximately constant total density (upper, blue 
curve). This behavior is reminiscent of dark soliton for- 
mation in colliding single-component BECs [18|. Beat- 
ing dark-dark solitons are also theoretically predicted to 
develop when a repulsive beam is swept through a two- 
component miscible BEC with an appropriate speed [19| . 

The beating solitons can be understood through the 
following simplified model: assuming that all scattering 
lengths are equal to <222, the mean field equation is the 
repulsive, vector Nonlinear Schrodinger (NLS) equation. 
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FIG. 5: (Color online) Dark-dark solitons as a result of two 
component mixing, (a) Phase separated mixture in presence 
of axial gradient, (b) Dark-dark soliton formed 1 sec after 
sudden turn-off of gradient, (c, d) Integrated cross sections 
with thin red (thick black) curve showing the |2, — 2) (|1, — 1)) 
component. Blue (upper) curve in (d) shows total density. 



Its most general known soliton solution is the six pa- 
rameter dark-dark soliton [4] (e.g. two background densi- 
ties 77,1,2, two background flow speeds ci ? 2, soliton speed 
v, and beating frequency uj) of which the well- studied 
five parameter static dark-dark soliton [3[ is a special 
case. Even though analytical expressions for these soli- 
tons were derived [4J , their form is quite complicated and 
basic properties such as the beating frequency as a func- 
tion of soliton parameters are unknown. 

An example of a beating dark-dark soliton can be con- 
structed by lev erag ing the SU(2) invariance of the vector 
NLS equation [15j. Applying a rotation matrix to the 
two components of a four parameter dark-bright soliton 
[3[, we obtain a five parameter beating dark-dark soli- 
ton where the background flow speeds are equal to c. Its 
evolution over half a beating period is shown in Fig. [6] 
(compare with Figs.[3fd,e) and[4je)). The beating angu- 
lar frequency uj = S(c — v) 2 sec 2 ((/>/2) satisfies [15| 

m(c — v) 2 /(2ft) < uj < 7rha22(ni + n^jm. (1) 



The soliton half- width is I = h/ \/2mujh — m 2 (c — v) 2 , 
where <j> is the soliton phase jump and m is the particle 
mass. As uj approaches the lower (upper) bound in ([2]), 
the beating soliton degenerates to a plane wave (static 
dark-dark soliton). The beating soliton strongly resem- 
bles features observed in experiment and numerical sim- 
ulations. The predicted minimum oscillation period of 
5 ms for our experimental parameters is consistent with 
the numerically observed periods of about 15 ms. 

In conclusion we have presented the first experimental 
observation of a beating dark-dark soliton. These solitons 
naturally arise from a fast counterflow-induced modula- 
tional instability and can emerge during the mixing of 
two superfluids. Our work opens the door to a range 
of new studies of vector soliton dynamics, with conse- 
quences for a diversity of nonlinear, dispersive systems. 
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FIG. 6: (Color online) Density and phase evolution of a beat- 
ing dark-dark soliton assuming equal scattering lengths. 
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SUPPORTING MATERIAL FOR "DARK-DARK SOLITONS AND MODULATIONAL INSTABILITY IN 

MISCIBLE, TWO-COMPONENT BOSE-EINSTEIN CONDENSATES" BY M. A. HOEFER, J. J. CHANG, 

C. HAMNER, AND P. ENGELS: NUMERICAL SIMULATIONS 

Numerical Computations 

Our 3D numerics are performed on the vector Gross-Pitaevskii equations 

ih-^- = "^ A *i + ^i( x ^)*i + (5n|*i| 2 +5i 2 |* 2 | 2 )*i, 

ih ^W = "£ A * 2 + ^ 2(X,f) * 2 + (^2l|*l| 2 +ff22|* 2 | 2 )*2, (2) 



f , T / mo , , T 4:7rh 2 aji 

/ |^(x,£)| 2 dx = 7V„ g 3l = ^, l<j,/<2. 

Jr3 m 



where gji = 47rh 2 aji/m, a,ji are the s-wave scattering lengths, m is the particle mass, Nj are the number of condensed 
atoms in each component, and 

^jCM) = 2 m ^ x2 +UJ v y2 +U3 lz 2 ) +VB9jmjB'(t)z, j = 1,2, (3) 

with the time varying magnetic field gradient B'(t), \±b is the Bohr magneton, the Lande g- factor is gj which has 
opposite signs for each component g\ = —1/2 and #2 = +1/2, hence field gradient induced counterflow. The hyperfine 
quantum number is rrtj. 

We use the split-step pseudospectral Fourier method [1] (adapted to multiple components) with a uniform grid 
spacing of approximately 0.19 jam in a rectangular 11.7 /xm x 11.7 /xm x 772 /xm box with periodic boundary conditions 
and a time step of approximately 0.0089 ms. The initial condition is computed by iteration of the stationary equations 
leading to the ground state in the presence of the time-independent portion of the trap. The experimentally invoked 
free expansion directly before imaging the condensate was not performed in the numerical simulations. As is common 
with the numerical simulation of unstable systems, we find that the onset time of modulational instability depends on 
the amount of noise in the system. Therefore, we construct noise with zero mean, 10 -4 variance, Gaussian distributed 
Fourier components for the smallest 32 wavenumbers along each spatial dimension. (1 + noise) multiplies the initially 
smooth density and is used as the initial condition. 

The movies contained in the EPAPS material show results obtained from our 3D numerics, simulating the ex- 
perimental procedure of Figs. 1 and 3 of the main text. The files contour _st rong_gr ad ient .large .mov and 
contour _strong_gradient_zoom_large .mov correspond to the experiments in Fig. l(a-c) and the numerical simu- 
lations in Fig. 4(b,d,f) of the main text. The file contour jnoder ate .gradient .large, mov corresponds to the ex- 
periments in Fig. 3(a,b,d,e) and the numerical simulations in Fig. 4(a,c,e) in the main text. In all cases, smooth 
counterflow leads to the spontaneous generation of a number of beating dark-dark solitons and their decay due to 
transverse modulations. 

Beating Dark-Dark Soliton 

The two-component vector Nonlinear Schrodinger equation modeling the (l+l)d dynamics of untrapped, binary 
BECs with equal scattering lengths is 

#t = -^„ + II^H 2 ^, il> = {tl>i,th) T - (4) 

This equation can be obtained from eq. © by a suitable dimensional reduction (see e.g. |2[) and nondimensionalizing 
lengths by the transverse harmonic oscillator length clq = ^/h/muj xi time by the transverse trap frequency lj x , and 



density by 2i\a s a^ where we have assumed that all scattering lengths are equal to a s . Four parameter dark-bright 
soliton solutions are well-known [3[. Because eq. (j4]) exhibits SU(2) invariance, we can "rotate" the most general 
dark-bright soliton and obtain a five parameter beating dark-dark soliton which was discussed in [4] 



^1 = yJ~P\ {cos(/> + isin(/>tanh[a(z — vt)}} ex.p \i[cz — (c 2 /2 + pi + p2)t]} — 

h[a(z - vt)} exp {i[vz + (a 2 /2 - v 2 /2 - pi - p 2 )t}} , 



P2SU1 



2 P 2 
a z secli 



Pi +P2 
^2 =\/~p2 {cos^ + isin^tanh[a(z — vt)]}exp {z[cz — (c 2 /2 + pi + P2)^]} + 



(5) 



pi sin 



2 P 1 
a z sec. 

pi + P2 



h[a(z - vt)] exp {i[vz + (a 2 /2 - v 2 /2 - pi - p 2 )t]} . 



The soliton parameters are the far field densities \^j\ 2 — » pj, |z| — >• 00, the soliton inverse half- width a, the phase 
jump 2(j) is the same for each component, the soliton speed v, the background flow speed c which is also the same for 
each component, and the beating frequency uo due to the relative, time dependent phases multiplying the tanh and 
sech terms. There are five independent parameters with the other two, say a and (/>, related through 



r- , a = \/2uj — (c — v) 2 . 



(6) 

Due to parameter restrictions on the original, un-rotated dark-bright soliton, the dark-dark soliton exists only for 
beating frequencies in the range 



(C-V) 2 gl+ p2 

<. U) <. . 



(7) 



One can directly verify that the beating dark-dark soliton in (J5j) bifurcates from a plane wave solution at w = (c — v) 2 /2 
and a static dark-dark soliton for uj = (pi + p2}/2. 
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FIG. 7: Sound speeds as functions of the relative speed. The smallest speed becomes purely imaginary for relative speeds above 
critical v cr , indicating the onset of modulational instability. 



For completeness, we include a plot (Fig. [7)) of the axial sound speeds [5[ 
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(8) 



ra z 4 x m* aum* 

for the case when the densities and scattering lengths satisfy the relation 

aii|^i| 2 = a 2 2|^2| 2 , (9) 

which corresponds to the experiments and analysis in Fig. 2 of the main manuscript. The relative superfluid speed is 

ft 



^rei = — |(arg* 2 )*-(arg*i)s|. 



m 



(10) 



Electronic address: mahoefer@ncsu.edu 

[1] W. Bao, D. Jaksch, and P. A. Markowich, J. Comp. Phys. 187, 318 (2003). 

[2] P. G. Kevrekidis, D. Frantzeskakis, and R. Carretero- Gonzalez, Emergent Nonlinear Phenomena in Bose-Einstein Conden- 
sates (Springer, Berlin, 2009). 

[3] A. P. Sheppard and Y. S. Kivshar, Phys. Rev. E 55, 4773 (1997). 

[4] Q.-Han Park and H. J. Shin, Phys. Rev. E 61, 3093 (2000). 

[5] C. K. Law, C. M. Chan, P. T. Leung, and M.-C. Chu, Phys. Rev. A 63, 063612 (2001). 



